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\ We numerically construct an one-parameter family of initial data of an expanding 

inhomogeneous universe model which is composed of regularly aligned black holes 
CN ■ with an identical mass. They are initial data for vacuum solutions of the Einstein 

\ equations. We call this universe model the "black hole universe" and analyze the 

structure of these initial data. We study the relation between the mean expansion 
rate of the 3-space, which corresponds to the Hubble parameter, and the mass density 

X- 

5_i , of black holes. The result implies that the same relation as that of the Einstein-de 

Sitter universe is realized in the limit of the large separation between neighboring 
black holes. The applicability of the cosmological Newtonian iV-body simulation to 
the dark matter composed of black holes is also discussed. The deviation of the 
spatial metric of the cosmological Newtonian iV-body system from that of the black 
hole universe is found to be smaller than about 1% in a region distant from the 
particles, if the separation length between neighboring particles is 20 times larger 
than their gravitational radius. By contrast, the deviation of the square of the 
Hubble parameter of the cosmological Newtonian iV-body system from that of the 
black hole universe is about 20% for the same separation length. 
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I. INTRODUCTION 

The homogeneous and isotropic universe model has enjoyed great success in explaining 
the observational data. By contrast, as anyone well knows, our universe is not exactly 
homogeneous and includes a lot of objects which serve as local nonlinear inhomogeneity. 
Usually, effects of local non-linear structures on the global property of the universe are 
considered by intuitive way or using some approximate methods. One of the effective ways 
to test the validity of our intuition or the approximation is to construct and study an exact 
or almost exact solution of the field equations, which may not be so realistic but should be 
able to fully describe non-linear effects in an inhomogeneous universe model. 

One example of exact inhomogeneous solutions is the so-called Swiss-cheese universe 
model [H, Q : the dust in arbitrary number of non-overlapping spherical regions is removed 
in a model of the homogeneous and isotropic universe filled with dust, and then each removed 
region is filled with a Schwarzschild black hole of the same mass as that of the removed dust. 
The remaining dust filled region, which is corresponding to "cheese" , is playing a role of the 
glue to connect Schwarzschild patches. However, due to the existence of the cheese region, 
the Swiss-cheese model may be too special to see significant effects of local inhomogeneities 
on the global evolution of the universe. Hence, it is important to study a universe model 
in which black holes are uniformly distributed without the cheese region. We call such an 
inhomogeneous universe model the "black hole universe" in this paper. 

About this issue, one innovative work has been done by Lindquist and Wheeler in 1957j3| 
and this work has been recently revisited in Refs. 0, \^ . They divided a virtual 3-sphere into 
N cells (N =5, 8, 16, 24, 120 and 600) and put a black hole portion of the Schwarzschild 
spacetime on a spherical region centered in each cell. Then they derived the equation 
of motion for this "lattice universe" from junction conditions between the Schwarzschild 
cell and the 3-sphere. It is demonstrated that the maximal radius of the lattice universe 
asymptotes to that of the corresponding homogeneous and isotropic closed universe filled 
with dust in the limit of the large number of black holes. Here we should note that the lattice 
universe is not an exact solution and there are gaps between each Schwarzschild black holes 
(see a figure 3 in Ref. j^]). 

Our purpose in this paper is to numerically construct initial data of the black hole uni- 
verse. As a first step, we consider regularly aligned black holes with an identical mass. By 
its symmetry, no anisotropic relative velocities between neighboring black holes will appear, 
and this system is similar to a cold gas, i.e., dust. In order to obtain such initial data sets, 
we consider a black hole at the center of a cubic region and impose the periodic boundary 
conditions on its faces. A recipe for the initial data of the black hole universe and numer- 
ical procedure are given in Sec. [Ill The degree of inhomogeneity of the black hole universe 
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is demonstrated in Sec. Ill I Al by calculating the traceless part of the 3- dimensional Ricci 
curvature tensor of the initial hypersurface. The structure of the initial hypersurface is 
investigated in Sec. IIIIBl by searching for horizons. 

One of the fascinating issues of inhomogeneous universe models is the so-called averaging 
problem. Naively, we expect that a universe model with local inhomogeneities, such as the 
black hole universe, can be globally described by a homogeneous universe model on average. 
However, the effect of local inhomogeneities on the global expansion definitely exists and the 
expansion history may be different from that of the homogeneous and isotropic universe 

0- 



9[. This issue has been discussed a lot in past years(see reviews [10Hl2[ and references 
therein), however there are few analyses which are applicable to inhomogeneous models with 
highly nonlinear metric inhomogeneity. 1 To solve this issue, we need to rely on numerical 
relativity. Our work may be the first step for the concrete analysis of the effects of non-linear 
inhomogeneities in expanding universes. Although we cannot address the real time evolution 
of the black hole universe yet, the one-parameter family of initial data sets can be regarded 
as a fictitious time evolution of the black hole universe. Using the initial data sets, we study 
the cosmic volume expansion rate of the black hole universe model in Sec. IIII Ci 

The cosmological A^-body simulation is a powerful tool for studying the structure forma- 
tion in the universe by dealing with the motion of point particles, based on the cosmological 
Newtonian approximation. Since the interaction between these particles is the only gravity, 
the cosmological A^-body simulation follows the time evolution of the dark matter in the cos- 
mological context. The black hole is a candidate of the ingredient for the dark matter, and 
it is believed that the cosmological A^-body simulation is applicable also to the black-hole 
dark matter. But it is quite non-trivial issue whether the point particles in the cosmological 
A^-body simulation can be simply identified with black holes. Hence, it is important to see 
in what situation the cosmological Newtonian A^-body simulation is valid for the black hole 
universe. This issue is discussed in Sec. IIII Dt Sec. [IV] is devoted to a summary. 

In this paper, we use the geometrized units in which the speed of light and Newton's 
gravitational constant are one, respectively. 



1 One of few exceptional examples was given in Ref. (l^ . They studied the volume expansion rate of a kind 
of the Swiss-cheese model and showed that the cosmic volume expansion can be accelerated by non-linear 



inhomogeneities. While we were writing this paper, Ref. [14[ appeared. In Ref. [14[, the authors analyt- 
ically constructed iV-body solutions of Einstein's constraint equations by considering regularly arranged 
distributions of discrete masses in topological 3-spheres. Significant differences between our present work 
and Ref. 1J| is the spatial topology and the existence of the cosmic volume expansion. In our present 
case, the spatial topology is T 3 with one point removed, and the expansion rate is finite while the initial 
date sets considered in Ref. Q have a topology of S 3 with N points removed, and their expansion rates 
vanish, i.e., time symmetric. One of remarkable advantages of our work over Ref. is that a dynamical 
simulation of expanding universe is possible starting from our initial data, while only contracting universe 



is possible with the initial data given in Ref. [14 [ 
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II. CONSTRUCTION OF INITIAL DATA FOR THE BLACK HOLE UNIVERSE 



A. Constraint equations 

In this paper, we are interested in the initial data of the vacuum Einstein equations. The 
initial data of the Einstein equations is equivalent to intrinsic and extrinsic geometries of a 
spacelike hypersurface, i.e., the intrinsic metric 7^ and the extrinsic curvature Kij, which 
represents how the spacelike hypersurface is embedded into the 4-dimensional spacetime. 
These are partially determined by the following four components of the Einstein equations: 



K + K 2 - K^KV = 0, 
DjK\ - DiK = 0, 



(1) 
(2) 



where 1Z and Di are Ricci curvature scalar and the covariant derivative with respect to the 
intrinsic metric 7^, respectively, and K = j^K^. Equation (JTJ) is called the Hamiltonian 
constraint, whereas Eq. (j5J) is called the momentum constraint. 



Following an established procedure (see, e.g., Ref. [l5|), we adopt the Cartesian spatial 
coordinate system and rewrite 7^ and as 



lij ^ Jij 



-10 



D l X 3 + D 3 X l 



f j D k X k + A ij 



TT 



(3) 
(4) 



where ^ := (det 7^)12, Di is covariant derivative with respect to the conformal metric 7^, 



and satisfies the transverse and traceless conditions, 



, jijA 



1,1 

TT 



0. 



(5) 



The conformal factor \1/ is determined so that the constraint equations are satisfied. The 
conformal metric 7^ has not six but five independent components due to the constraint 
det7ij = 1. The three of the five components of 7^ can be always eliminated by the spatial 
coordinate transformation, and hence there are two physically meaningful components which 
can be freely chosen. 

In the decomposition (j4j), mutually independent six components of are expressed by 
X 1 , A^ T and K. The longitudinal traceless part composed of X 1 is determined so that 
the constraint equations are satisfied, whereas the trace part K is related to the degree of 
freedom to choose the foliation of the spacetime by the family of spacelike hyper surf aces, 
or in other words, time slicing. By contrast, the transverse and traceless part has two 
independent components which can be freely chosen. These two components of and 
the physically meaningful two components of 7^ are usually regarded as physical degrees of 
freedom to set initial data for gravitational waves. 

In order to avoid the cosmic volume expansion caused by artificial gravitational radiation, 
we assume trivial form of the conformal metric and no transverse and traceless part of the 
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extrinsic curvature 

lij = Sij, (6) 
4 T = 0, (7) 

where 8^ is Kronecker's delta. As usual, we denote the inverse of 7^ by 7* J which is also 
equal to Kronecker's delta 8 l i . Then, Eqs.Q and (j2j) are written as 

A* + l(LX) tJ (LX) ij V 7 - ^-K 2 ^ 5 = 0, (8) 
o 12 

AX 1 + -&djX j - = 0, (9) 

3 3 

where A is the flat Laplacian, di is the ordinary derivative, and 

(LX) ij := d i X j + d j X i - -5 ij d k X k . (10) 

3 

Here, note that X % := j^X* = X i and & := Y j dj = di, (LX)^ = ^ ik j jl (LX) kl = (LX)*, 
etc. We solve these equations by assuming an appropriate functional form of K as shown 
below. 



B. Boundary condition and the trace of the extrinsic curvature 

As mentioned above, we adopt the Cartesian coordinate system x = (x, y, z) and put 
a non-rotating black hole at the origin x = denoted hereafter by O. The black hole is 
represented by a structure like the Einstein-Rosen bridge in our initial hypersurface. Thus 
the origin O corresponds to the asymptotically flat spatial infinity and is often called the 
puncture. We focus on a cubic region — L < x < +L, —L < y < +L and — L < z < +L and 
call this region the domain T>. Since our aim is to construct the initial data of an expanding 
universe model with periodically aligned black holes, we impose the periodic boundary 
conditions; a point x = (—L,y,z) is identified with a point x = (+L,y,z), etc. Due to 
this boundary condition, the domain T> is homeomorphic to the 3-torus T 3 . Since infinity 
is not included in the spacetime manifold, the initial hypersurface is T> with O removed, 
which is denoted by D — {O}, and thus it is homeomorphic to T 3 with one point removed. 2 
The covering space of T> — {0} represents a cosmological model with regularly aligned black 
holes as shown in FigJTJ Hereafter, we regard P as a cubic domain with boundary &D in 
the covering space. 

Here, we again note that the trace part of the extrinsic curvature K corresponds to the 
degree of freedom to choose the time slicing. In order to find the appropriate time slicing 
condition, first of all, we see the homogeneous and isotropic universe model. In this case, the 



2 A similar configuration to our case was considered within the Lemaitre-Tolman family of exact models 
in Refs. 0, 
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FIG. 1: The cubic region of our coordinates. 



expansion rate H which is called the Hubble parameter is related to the extrinsic curvature 
by 

H=-~K. (11) 
3 



The above relation implies that K of the expanding black hole universe model must be 
negative at least around the boundary of the cubic domain T>. By contrast, the maximal 
slicing condition K = is appropriate for the foliation of the domain in the neighborhood 
of the asymptotically flat spatial infinity, and hence K should vanish in the vicinity of O 
(see AppendixTAl). 

Taking the above discussions into account, we assume 



K(x] 



-3H eS W(R), 



(12) 



where if e ff is a positive constant which corresponds to the effective Hubble parameter, R 
Ice I, and 



W(R) 



-36 



[(R 



a 



for < R < i 
a e ] 6 for t < R < i + a 
for i + a < R 



(13) 



and a being constants which satisfy £ < a < L (see FigfS]). 
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FIG. 2: The functional form of W(R). 



C. Extraction of the singularity at the center 



Since K vanishes in the vicinity of the origin O, X % and ^ should behave as those of the 
Schwarzschild spacetime with the static isotropic coordinate system, 



X 1 



-o, 

M 



(14) 
(15) 



where \l/ c and M are positive constants. Since \1/ itself is singular at O, we cannot handle 
\1/ numerically. Thus, instead of we solve the constraint equations for the following new 
variable ip: 

M 

(16) 



rf>(x):=V(x)- — [l-W{R)]. 

Thanks to the second term proportional to [1 — W(J?)] in the right hand side of the above 
equation, ip is regular at O and satisfies the periodic boundary conditions. 

The mass of a black hole is given by the ADM mass which is defined by the surface 
integral over the spacelike infinity at O. To see the ADM mass explicitly, we introduce a 
new radial coordinate 

~ M 2 

R -lR- < 17 > 
Then, by using a spherical polar coordinate system, the asymptotic form of the infinitesimal 
line element for R — y 0, or equivalently, R — » oo, becomes 



v& c + ||) [di? 2 + i? 2 (dfl 2 + sin 2 #d0 2 )] 



1 + 



2R 



dR 2 + R 2 (d9 2 + sin 2 



(18) 
(19) 



s 



It is seen from the last equality of Eq. (I19p that the mass of a black hole is given by \I/ C M. 
Here note that there is a freedom of constant scaling of coordinates x — > Cx. Using this 
freedom, we impose \l/ c = 1, and thus the mass of a black hole is equal to M. 



D. Hubble equation as an integrability condition 

Integrating Eq.flEl) over the physical domain T> — {O}, we obtain the following equation: 

2ttM + I [ {LX^LXf^dx 3 - -H 2 cS V = 0, (20) 
S Jv-{o} 4 

where, by noting that the origin O can be regarded as the only boundary of D — {0} with 
the periodic boundary condition, the integral of A\l/ is rewritten as 

r r 2lT r 71 f)\tj 

/ A*d 3 x = -lim/ / — R 2 sm9d$d(f) = 2irM, (21) 
Jv-{o} ™Jo Jo dR * y J 

and we have defined V by 

V := [ W 2 ^ 5 d 3 x. (22) 
By rewriting Eq. (120]) . we have the effective Hubble equation as 



57T 



#cff = y (Pbh + Pk) , (23) 



where pbh and pk are defined by 



M , . 

Pbh := y, (24) 

Pk := -^TT / (LX) tJ (LXy^- 7 d 3 x. (25) 



16nV 



v-{o} 



Since V may be regarded as the effective volume of the expanding region, p B H and px m ay be 
regarded as the mass density of black holes and the kinetic energy density of the spacetime, 
respectively. The effective Hubble equation gives a relation between two constants, the mass 
of the black hole M and the effective Hubble parameter H e g, and, at the same time, it is an 
integrability condition of the constraint equations. How to guarantee this relation will be 
described in Sec. Ill Fl 



E. Momentum constraints 



In this subsection, we rewrite the momentum constraints (Q into the numerically solvable 
forms. First, we define Z by 

Z := diX 1 . (26) 
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Then, by taking the divergence of Eq. (Q, we obtain 

AZ = -diiVt&K). (27) 

Eq. (Q is rewritten as 

AX 1 = -d l Z + -V*&K. (28) 
3 3 

The system we consider is unchanged if it rotates 27r/3 radians around the line x = y = z. 
By virtue of this discrete symmetry, it is enough to solve Eq. f[2"8"j) for only one component 
of X 1 , since the other two components can be immediately given by this symmetry. 
The boundary condition for X x is given as follows, 

X x = on x = and x = L , (29) 
d y X x = on y = and y = L , (30) 
d z X x = on z = and z — L. (31) 

The first condition is the Dirichlet type, and the second and third ones are Neumann type 
boundary conditions. These boundary conditions lead to 

/ Zdx 3 = [ d i X i dtf s = 0. (32) 
Jv~{o} JV~{0} 

The above equation is a consistency condition that the solutions should satisfy. 

It should be noted that the integrals of the source terms of the Poisson equations (I27|) 
and (128]) over D — {0} should vanish by the consistency with the periodic boundary con- 
ditions and the boundary condition at the origin O. We can see that these conditions are 
automatically satisfied. The integral of the source term of Eq. (1271) is equivalent to the 
surface integral over the spatial infinity at O, whereas K vanishes in the neighborhood of 
O. Hence the integral of the source term of Eq. (J27J) vanishes. Since d x Z and ^> 6 d x K are 
odd functions of x, we have 

" +L f 4 1 \ 

dx f ~^d x Z + -^ 6 d x K J = 0. (33) 

Hence, the integral of the x-component of the source term of Eq. (128!) vanishes. The same 
is true for the other components of Eq. (128]) . 



F. Numerical procedure 

As shown in the preceding section, we have to solve the following three coupled Poisson 
equations, 

= A (^W(R)^ -±{LXMLX)*V- r + ±K*V*, 

AZ = ^(^W), 
AX 1 = --&Z + \^d l K. 
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In order to get numerical solutions of the above equations, we adopt the method of fi- 
nite differentiations. By replacing all derivative terms by finite differences, we have a very 
large simultaneous equation. We solve this simultaneous equation by the Successive Over- 
Relaxation method. We denote the values of x/i, Z and X % at each iteration step by ipQ, 
ifti, ip2 and so on, where the subscript denotes a trial value. At the (n + l)-th step 
of the iteration, the terms corresponding to the source terms of the Poisson equations are 
estimated by using ip n , Z n and X l n . 

If we complete the n-th step of the iteration, we obtain Z n which satisfies the boundary 
conditions (l29l)-( |3T|) . Here, we should note that this Z n does not necessarily satisfy the 
consistency condition ( 1321) . In order to obtain Z n which satisfies Eq. (I32jl . we can use the 
degree of freedom to add a constant to Z as follows, 

Z ->• Z' := Z - / dx 3 Z. (34) 

Z' is also a solution of Eq. ( |27l) and further satisfies Eq. ( )32|) . if Z is a solution of Eq. ( |27l) . 
Thus, before evaluating the source term, we reset the value of Z n as follows: 

Z n ^Z' n = Z n -±- [ dx 3 Z n . (35) 
L Jv-{o} 

It should also be noted that the boundary conditions already given are not enough to close 
the simultaneous equation, since these boundary conditions do not determine homogeneous 
solutions of the Poisson equations for ip and X 1 , i.e., their zero modes. (The zero mode of 
Z is already fixed by Eq. (|35|) .) For this purpose, we need to specify the values of if) and 
X 1 at one of all numerical grids. We fix the zero modes of ip and X % so that ip(0) = 1 and 
X*(0) = 0, or in other words, before evaluating the source terms, we add constants to ip n 
and X l n as 

i> n (x) — ► tf n {x) := 4> n {x) - ^ n (0) + 1, (36) 
X* n (x) — ► X'iix) := Xl{x) - X;(0). (37) 

Note that ^(0) = 1 is equivalent to the choice of \l/ c = 1 in Eq. ( fl5l) . Eventually, we evaluate 
the source terms by using ip' n , Z' n and X' l n instead of ip n , Z n and X % n . The value of if e ff is 
also evaluated through Eq. ( f23l) by using tf}' n , Z' n and X'\ so that the integrability condition 
is satisfied. 



G. Results 

We solved the constraint equations in the parameter domain 2.8i? g < L < 20R g , where 

R g := f. (38) 

As will be shown later, the horizons of a black hole are located at R ~ R g . The parameters 
a and t which determine K are set to be a = 0.2R„ and £ = L — 0AR„. We could not 
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get convergence for L smaller than 2.8R g . This result implies that there is no solution for 
L < 2.8R g on our assumptions: conformally flat metric and no transverse-traceless part of 
the extrinsic curvature. The results of convergence test for each value of L/R g is shown in 
Fig. |3j The second order convergence is confirmed in all cases for the value of where 
the reference value i? r 2 cf is given by the least-square fit. We plot ip, Z and X x on z = and 
z = L planes as functions of x and y for L — 2.8R g and L = 10R g in Figs. HI [5] and |6j 



3) 

CM 




5e-005 



0.0001 

(dx / L) 2 



0.00015 



0.0002 



FIG. 3: Results of convergence test. 



<|/ on z=0 for L=2.8R, 




y on z=L for L=2.8R, 



1.44 r 




>|/on z=Lfor L=10R, 



1.135 r 




FIG. 4: tp on z = and z = L planes as functions of x and y for L = 2.8R g and L = lOii, 
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Z[R ]on z=0forl_=2.8R 



Z[R g - 1 lonz=LforL=2.8R g 





Z[R "']on z=0 for L=10R 



Z[R "']on z=Lfor L=10R 





FIG. 5: Z on z = and z = L planes as functions of x and y for L = 2.8R g and L = 10R„. 



X on z=0 for L=2.8R. 



X" on z=L for L=2.8R„ 
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X x on z=0forL=10R. 



X x on z=LtorL=10R„ 





FIG. 6: X x on z = and z = L planes as functions of x and y for L = 2.8R g and L = 10R„ 
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III. ANALYSIS OF THE INITIAL DATA 



A. Inhomogeneities 

First, we demonstrate the inhomogeneities of our initial data. For this purpose, we 
investigate the following quantity: 

_ T c i bd K b nl 

where IZij and TZjj denote the 3-dimensional Ricci curvature tensor and its traceless part, 
respectively. We use (3 as a measure of homogeneity and isotropy, since a region with (3 = 
and d{IZ = is homogeneous and isotropic. Since we are interested in the inhomogeneities 
far from black holes, we plot the value of on z = L plane, which is one of the faces 
of the domain V, as a function of x and y in Fig. [7J The quantity /3 almost vanishes in 
the vicinity of a vertex x = y = z = L. Further, the norm of the traceless part of the 
extrinsic curvature fy~ l2 (LX)ij(LXy : > is much less than the square of the trace part of the 
extrinsic curvature K 2 in the neighborhoods of the vertices. We find from the Hamiltonian 
constraint together with this fact that 1Z ~ —K 2 =constant and hence IZij ~ —3H 2 e 'jij, in 
these regions. Thus, the neighborhoods of the vertices are well approximated by the Milne 
universe model which is the Minkowski spacetime foliated by the family of homogeneous and 
isotropic spacelike hypersurfaces with negative Ricci curvature scalar. Conversely, around 
the center of a face of T> {x = y = and z — L), the inhomogeneity remains even if L 3> -R g . 
We may understand this result as follows. If the neighborhoods of all faces of T> were well 
approximated by the Milne universe model, a 3-hyperboloid would be tiled with the lattice 
structure shown in Fig. [TJ However, this consequence conflicts with a mathematical theorem 
about "tiling" 0, 18]. Therefore, our initial data cannot be homogeneous and isotropic in 



the neighborhoods of all faces of T> even for L ^> R„ as is explicitly shown in Fig. [7J 



B. Horizons 

We define a horizon as a spacelike closed 2-surface with vanishing expansion of a null 
vector field normal to the 2-surface. There are two independent null directions normal to 
the 2-surface, so there are two kinds of horizons accordingly. Here, we consider these horizons 
in the domain T> — {O}. A closed 2-surface divides the domain V> — {0} into two regions. In 
this paper, since we are interested in the horizons associated with a black hole, we focus on 
a case in which one of the two regions includes the puncture. We call the domain including 
the puncture the inside, whereas the other domain is called the outside. Then, we call the 
direction from a point on a closed 2-surface to the outside the outgoing direction, whereas 
the opposite direction is called ingoing direction. Accordingly, we name a horizon with 
vanishing expansion of the outgoing null vector field the black hole (BH) horizon, whereas 
a horizon with vanishing expansion of the ingoing null vector field is named the white hole 
(WH) horizon. 
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The expansions of the null vector fields normal to this 2-surface are given by 

X± = (7 ij '-«V)(±Asi-#y), (40) 

where the subscript + means that of the outgoing null, whereas the subscript — represents 
that of the ingoing null, and s l is the outgoing unit vector which is normal to this 2-surface 
and tangent to the initial hypersurface. Defining s l and Sj as 



we rewrite x± i n t ne form 

X± = (siSj - 5ij) 



m~\LX) ij + -S ij K 
3 



(41) 



(42) 



In this paper, in stead of solving the equation x± — 0, we investigate the expansions of 
the null vector fields normal to various spheres centered at the origin O. The conformal unit 
vector s l normal to the sphere of the radius R is given by 



(43) 



If the initial hypersurface is almost spherically symmetric near the horizon, the horizon is 
also almost spherically symmetric and s l is a good approximation of the unit vector field 
normal to the horizon. In Fig. [HI we plot the expansions x± as functions of R on the following 
three lines: 



w y 

(ii) x 

(iii) x 



0, z = 0, 

y, 2 = 0, 

y = z. 



15 



From these figures, we see that there are spheres which are very good approximations 
of horizons; the expansion x+ 01 X- a t the intersections of these spheres and the lines (i) 
- (hi) vanishes. The coordinate radius R of the BH horizon is equal to 1.14R g in the case 
of L = 2.8-Rg, whereas it is equal to R g in the case of L = 10R g . The coordinate radius R 
of the WH horizon is equal to 0.92i? g in the case of L = 2.8R g , whereas it is equal to R g 
in the case of L = 10R g . In the case of L = 2.8R g , the WH horizon is located inside the 
BH horizon. Since the domain R < R g is well approximated by the Schwarzschild BH, we 
can say that the initial hypersurface is passing through the future of the bifurcation point 
of the horizons for L = 2.8R g . On the other hand, in the larger L cases, since the WH and 
BH horizons coincide with each other, we may say that the initial hypersurface is passing 
through a domain very close to the bifurcation point. 

XforL=2.8R XforL=10R. 
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FIG. 8: x± as functions of R for L = 2.8R g and L = WR i 



C. Effective Hubble equation 

The mass density of black holes pu defined by Eq. ([24]) is roughly estimated at about 
M/8L 3 . If the kinetic energy density px defined by ([25]) is much less than Pbh ; the effective 
Hubble parameter H c q is roughly estimated at about H^ s ~ 87t/7bh/3 ~ ttM/3L 3 . Then, in 
the covering space of the domain T> — {O}, the number TVbh of black holes within a sphere 
of the cosmological horizon radius H~g is about 

1 47r rr 3 1 ( 3 ^ 3 

N ™~M X T H ^~4{^J ' (44) 

If L/Rg is much larger than unity, there are many black holes within a sphere of the cos- 
mological horizon radius, and thus the black hole universe would be very similar to the 
Einstein-de Sitter(EdS) universe. 

From the above consideration, we expect that the effective Hubble parameter and the 
mass density of black holes asymptotically satisfy the Hubble equation of the EdS universe 
in the limit of L/R g — > oo. That is, we expect that the effective Hubble parameter behaves 
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asymptotically as 



8tt 



eff 



-PBH- 



(45) 



This means that the contribution of /?k decreases with larger L/R g . We depict Pk/pbh as a 
function of L/R g in Fig. |9j It is seen from this figure that Pk/pbh asymptotically vanishes 
for large L/R g and the effective Hubble equation approaches that of EdS universe. 
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FIG. 9: Pk/pbh as a function of L/R g . 



It is suggestive to regard the one parameter family of the initial data sets as a fictitious 
time evolution of the black hole universe. Eq. (|23|) gives the effective Hubble parameter at 
each time of the fictitious evolution. If we define an effective scale factor by 



ay 



V 1/3 , (46) 



Eq. (|45p means that H^ s asymptotically behaves as oc l/a v when the universe expands 
enough. 

Other remarkable ways to define effective scale factors are to use the proper area of the 
boundary and the proper length of the edge of the cubic domain T>. Let a a and ax denote 
the effective scale factors defined by using the proper area and the edge length, respectively. 

is defined by the proper length of a edge itself and a a is defined by 




aA:=\hr, (47) 



where A is the proper area of &D. In addition, we define the fiducial scale factor a^ds by 
using the Hubble equation of EdS universe as follows: 

4dS : = 17T7T- (48) 

^eff 
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The relation between effective scale factors and the effective Hubble parameter is shown in 
Fig. [lOj All effective scale factors asymptotically behave as oc H c ^ 3 for larger L/R g , that 



a L 3 /a Eds 3 -1 _ 

a A ' a EdS "1 ■"*■■■ _ 

3 3 
a V ' a EdS 



B 

» - 

., . . . 

100 1000 
1/(H e /R g 2 ) 

FIG. 10: Effective scale factors(left panel) and deviations of them from the fiducial scale factor 
o-EdS (right panel) as functions of the effective Hubble parameter. 

is, the behaviour of the effective Hubble parameter as a function of an effective scale factor 
agrees with that of the EdS universe at late time of the fictitious time evolution. We note 

2/3 

that, even though all effective scale factors are asymptotically proportional to H eS , the 
proportionality coefficients are different from each other. It seems that the proportionality 
coefficient for ay asymptotically agrees with that for a^dSi but h is n °t true for a a and 
at, (see the right panel of Fig. ITU]) . 

D. Comparison with the Newtonian approximation and backreaction effect 

One possible way of approximation which considerably reduces the numerical effort is the 
cosmological Newtonian approximation. The cosmological iV-body simulation based on this 
approximation scheme is very useful to study the structure formation in the universe indeed. 
The iV-body simulation follows the motion of point particles gravitationally interacting with 
each other, and these particles are regarded as the dark matter in the cosmological context. 
The black hole is a candidate for the ingredient of the dark matter in our universe. However, 
since the black hole is a highly relativistic object, it is non-trivial whether the dark matter 
composed of black holes is well described by the cosmological iV-body simulation based on 
the Newtonian approximation. Our black hole universe model is a relativistic version of the 
cosmological iV-body system, and thus, by using this model, we can see in what situation the 
Newtonian iV-body simulation correctly describes the motion of the dark matter composed 
of black holes. 

In the cosmological Newtonian approximation scheme, the gravitational force is given 
by the spatial gradient of the Newtonian potential $ which is related to the conformal 
factor \1/ by 1 — 2$ = \l/ 4 , and thus the Newtonian potential is obtained by solving the 
Hamiltonian constraint, which gives the Hubble equation after averaging. Since the metric 
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is assumed to be almost equal to that of the EdS universe model, the term proportional 
to {LX)ij{LX) 1 ^ should be so small that it is a negligible higher order correction in the 
Hamiltonian constraint. Hence, we do not need to solve the momentum constraint. 

Before considering a point particle as the ingredient of the iV-body simulation, we assume 
that the particle is a spherical ball with the finite energy density p(x). Further, we assume 
the similar situation to our black hole universe; the particle has the mass M, the center of 
the particle is located at the origin O in the cubic domain T> whose edge length is 2L, and 
the periodic boundary condition is imposed. By definition, we have 

M= I p(x)d 3 x. (49) 
Jv 

The time slicing condition up to the Newtonian order is assumed to be 

K = -3# N , (50) 

where is the effective Hubble parameter up to the Newtonian order and is determined 
by 

if 2 - ^ x 

3 8L3- (51j 
Here note that is the same as the Hubble parameter of the background EdS universe 
model. Then, since nonlinear terms with respect to \I> in the Hamiltonian constraint is 
linearized with respect to $, the Hamiltonian constraint takes the following form in the 



cosmological Newtonian scheme 19j, |20 



' M 



(52) 



In the cosmological Newtonian approximation scheme, p can be much larger than M/8L 3 , 
but p should be so small that |<fr| is much smaller than unity. 

Let us consider the case in which the size of the particle is much smaller than L. In 
this case, since the tidal force can be neglected, it is enough to consider the energy density 
for a point particle given by M5(x) instead of the finite energy density p{x). Using this 
approximation, we can accurately estimate the gravitational force produced by a particle 
at points of other particles. Then the Hamiltonian constraint in the cosmological iV-body 
system is given by 

7rM 

A$ = 4nM5(x) — . (53) 

ZIj 

Equation ( 153]) is the basic equation for the cosmological iV-body simulation based on the 
Newtonian approximation scheme. 

In our case, since the ^> diverges at O in the black hole universe, it is obvious that 
the cosmological Newtonian approximation is not applicable in whole region of T> — {O}. 
When we compare the black hole universe with the cosmological Newtonian system given 
by Eq. (1531 . the point-particle approximation, i.e., p(x) = M6(x), should be regarded as a 
technical simplification. Hence, it is a very non-trivial issue whether a point particle in the 
cosmological Newtonian iV-body system may be identified with a black hole. 
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In order to numerically obtain solutions of Eq. ( |53|) . we decompose $ as follows, 



^ = ( f ) -^[l-W(R)}. 



The equation for <fi is given by 



A0 = -A (f ^(i?)) 



irM 

21?' 



(54) 



(55) 



This can be numerically integrated in T> — {0} by using the same way as in Sec. Ill Fl 

To show the deviation of the solution obtained by the cosmological Newtonian approxi- 
mation from the corresponding relativistic one, we plot the following quantity: 



a 



2$ 



\]/4 



(56) 



In Fig. HU a is plotted as a function of the coordinates x and y on z = and z = L planes 
for L = 2.8i?g and L = 20R g , respectively. We can see that while the deviation around the 
boundary of T> is a few tens of percent for L = 2.8R g , it is less than 1% for L = 20R g . This 
result implies that the cosmological Newtonian approximation predicts the spatial metric 
around the boundary of T> very accurately for L > 20R g . 
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FIG. 11: a on z = and z = L surfaces for L = 2.8R P and L = 20R 
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As mentioned, the effective Hubble parameter defined by Eq. (ISTil agrees with that 
of the background Einstein-de Sitter universe model. The so-called backreaction effect is the 
change of the Hubble parameter from the background value due to the nonlinear effect of the 
inhomogeneities. Thus, in the present case, we call the effect which causes difference between 
the full relativistic Hubble parameter H eS and the background value H N the backreaction 
effect. 
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To see the significance of the backreaction effect, we compare H^ s to with fixed L/R s . 
As a result of the numerical investigation, we find that has about 20% deviation from 
ifg ff even for L = 20R g . We plot the value of l—H^/H^ as a function of L/R s in Fig. [121 It 
is worthwhile to notice that the Newtonian Hubble parameter is larger than the relativistic 
one. This means that the backreaction effect acts as the brake in the black hole universe 
model. Further, our result means that, in the case of L < 20R g , the backreaction effect is 
so large that the cosmological Newtonian approximation cannot predict correctly the global 
cosmic volume expansion rate. However, Fig. [12] suggests that the deviation of from 
if e ff decreases with larger L/Rg, and hence it seems that the Newtonian A-body simulation 
becomes correct asymptotically for L/Rg — > oo. 

As already shown, in the case of L = 20R g , the relative difference in the spatial metric 
between the Newtonian scheme and the full relativistic one is a few percents on the boundary 
of V, and hence the relative differences in the length of an edge and the area of a face are 
also a few percents. Furthermore, px defined by Eq. (J25l) is about 2% of pbh (see Fig. [9]). 
Thus, the difference between and comes from the difference between the volume V 
defined by Eq. (j22l) and 8L 3 ; V is about 1.3 times larger than 8L 3 . 

Here, we should note that the backreaction effect is large even in the case of L = 20i? g , 
but, as shown in the preceding section, the expansion law of the black hole universe model 
might be almost the same as that of the EdS universe. These results would imply that the 
backreaction effect in the black hole universe model would not change the expansion law 
from the EdS universe model but apparently shifts the time to the future. However, in order 
to get definite conclusion, the investigation of the time evolution is necessary. 
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FIG. 12: 1 - H^/H^ as a function of LjR x 



IV. SUMMARY AND CONCLUSION 

In this paper, we have constructed numerically the initial data of an expanding universe 
model which is composed of regularly aligned black holes. This system is equivalent to a 
black hole located at the center O of a cubic domain T> with periodic boundary conditions. 
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The black hole is represented by a structure like the Einstein-Rosen bridge, and thus O 
corresponds to the asymptotically flat spatial infinity. Since the physical domain does not 
include infinity, the physical domain is T> with O removed, i.e., T> — {0} whose topology is 
T 3 with one point removed. The functional form of the trace of the extrinsic curvature K(x) 
has been chosen so that K is a negative constant denoted by — 3H c g in the neighborhoods 
of the faces of T> and vanishes in the neighborhood of O, where H e g corresponds to the 
effective Hubble parameter. These requirements are compatible with a finite expansion rate 
of the universe and the puncture method to numerically treat a black hole, respectively. 
Then, we can solve constraint equations by giving the parameter L/R g , where L is the 
coordinate length of an edge of the cubic domain V, and R g gives a coordinate value which 
is almost equal to the coordinate radius of the black hole horizon. The value of H c g is 
determined so that the integral of the Hamiltonian constraint over T> — {0} is compatible 
with the periodic boundary conditions; this integral leads to the effective Hubble equation. 
We find from numerically obtained solutions that the neighborhoods of vertices of T> are well 
approximated by the Milne universe, whereas the other region remains inhomogeneous even 
in the case of L 3> -R g . This result implies that the initial data of the black hole universe 
model is inhomogeneous even near the faces of T> irrespective of the value of L/R g . 

We could find one white hole and one black hole horizons in the present initial hyper- 
surface of T> — {O}, and both are almost spherically symmetric. This result implies that 
the region R < R g is well approximated by the Schwarzschild black hole, and the initial 
hypersurfaces considered here are passing through the future of the bifurcation point of the 
horizons or a very close point to it. 

In order to compare our initial data with the Einstein-de Sitter(EdS) universe, we studied 
the relation between the effective mass density pbh of black holes and the effective Hubble 
parameter H e g, which are defined in a simple and natural way. Then, our numerical solutions 
imply that pbh and H eS asymptotically satisfy the Hubble equation of the EdS universe for 
L ^> R g . Once we regard our one parameter family of initial data sets as fictitious time 
evolution of the black hole universe, our result would imply that the Hubble equation of the 
EdS universe would be realized when the universe expands enough. 

The validity of the Newtonian approximation in the system has also been discussed. 
We numerically solved the Hamiltonian constraint equation simplified by the cosmological 
Newtonian approximation and compared it with the full solution with fixed L/R g . We 
found that the deviation of the spatial metric obtained by the cosmological Newtonian 
approximation from that of the full calculation is less than 1% for L/R g = 20 around the 
boundary of T> and better for larger values of L/R g . However, the deviation of the Hubble 
parameter defined in the cosmological Newtonian approximation scheme and full relativistic 
one is 20% even for L/R g = 20. Thus, we may say that, as expected, the backreaction 
effects of the inhomogeneities on the cosmic volume expansion is very large in the case of 
L < 20R g . However, we may also say that, for the larger L/Rg, the backreaction effects 
become smaller. It is worthwhile to notice that the backreaction effect acts as the brake for 
the cosmic volume expansion. 

Although our results in this paper agree with naive expectations, it is not clear by the 
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present analysis if the dynamics of the black hole universe can be described by the EdS uni- 
verse on average or not. Because the black hole universe cannot be exactly the EdS universe 
and the effect of inhomogeneities definitely exists. For instance, if we have an effective pos- 
itive curvature term on average as the effect of the inhomogeneities, the black hole universe 
eventually re-collapses. The effect of the inhomogeneities might give a qualitative difference 
of the global expansion history of the universe[8|, |9|, |2l|, |22j. By contrast, the present results 



would imply that the backreaction effect would not change the expansion law of the black 
hole universe from that of the EdS universe model; the backreaction effects might merely 
shift the time to the future. To attack this issue we need further numerical efforts and we 
leave it as a future work. 
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Appendix A: Constant Mean Curvature Slices in Schwarzschild Spacetime 

Let us consider the Schwarzschild spacetime, whose metric is given by 

ds 2 = -/(r)dt 2 + -^-dr 2 + rW, (Al) 
f{r) 

where 

f(r) = 1 - !l. (A2) 
r 

We consider a constant mean curvature(CMC) slice given by the form of 

t = h(r). (A3) 

The unit normal vector is given by 

= y, / fh', 0, )• (A4) 



Vf- 1 - f h ' 2 

The CMC slice condition is given by 
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45 —d r (r n 

1 C 



45 n 

^f- L (l-fh 



3 

f2u /2\ 



F(r;r g ,K,C) 



l-^ + (-iAV + g) 2 ' 

r y 3 r z ' 

where C is the integration constant. Then, line elements on the CMC slice is given by 

df = F(r; r g , K, C)dr 2 + r W. 
The transformation to the isotropic coordinate can be done as follows: 
df = f 4 (di? 2 + i?W), 



R = Cexp ± j dr^F(r;r g ,K,C)/r 

J r min 



(A5) 

(A6) 

(AT) 
(A8) 

(A9) 



where r min is the largest root of F(r; r g , K, C) = and r = r min corresponds to the throat. 
The minus sign is used in the region beyond the throat. We can easily check that, in the 
limit of r — > oo, the isotropic coordinate R has finite value if K ^ 0. While R = for 
r — > oo if K = 0. Hence, the coordinate region with R has inside spherical boundary with 
K ^ 0. This property is not compatible with the puncture method. 
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